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We perform numerical simulations to determine the shear stress and pressure of steady-state 
shear flow in a soft-disk model in two dimensions at zero temperature in the vicinity of the jamming 
transition <f>j. We use critical point scaling analyses to determine the critical behavior at jamming, 
and we find that it is crucial to include corrections to scaling for a reliable analysis. We find that the 
relative size of these corrections are much smaller for pressure than for shear stress. We furthermore 
find a superlinear behavior for pressure and shear stress above cj)j, both from the scaling analysis 
and from a direct analysis of pressure data extrapolated to the limit of vanishing shear rate. 
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Granular materials, supercooled liquids, and foams are 
examples of systems that may undergo a transition from 
a liquid-like to an amorphous solid state as some control 
parameter is varied. It has been hypothesised that the 
transitions in these strikingly different systems are con- 
trolled by the same mechanism [1] and the term jamming 
has been coined for this transition. 

Much work on jamming has focused on a particularly 
simple model, consisting of frictionless spherical particles 
with repulsive contact interactions at zero temperature 
[2]. The packing fraction (density) of particles </> is the 
key control parameter in such systems. Jamming upon 
compression, and jamming by relaxation from initially 
random states, have been the focus of many investiga- 
tions [2, 9, 11]. Another, physically realizable and impor- 
tant case is jamming upon shear deformation. This has 
been modeled both by simulations at a finite constant 
shear strain rate 7 [3-8, 10], as well as by quasistatic 
shearing [11, 15, 17], in which the system relaxes to its 
energy minimum after each finite small strain increment. 

Several attempts have been made to determine the crit- 
ical packing fraction </>,/ and critical exponents, describing 
behavior at shear driven jamming [3-8]. There is how- 
ever little agreement on the values of the exponents and 
there is thus a need for a thorough and careful investi- 
gation of the jamming transition in the shearing ensem- 
ble. It will also be interesting to compare the exponents 
found from shearing rheology to those found from com- 
pressing marginally jammed packings. In particular we 
note the linear increase of pressure above jamming that 
is observed in that system [2, 9], compared to the super- 
linear behavior often reported in the sheared system for 
pressure and/or shear stress [3, 4, 7, 8]. 

In this Letter we do a careful scaling analysis of high 
precision data for both shear stress and pressure at shear 
strain rates down to 7 = 10~ 8 . Instead of relying on 
visually acceptable data collapses we use a non-linear 
minimization technique to determine the best fitting pa- 
rameters. As in a recent analysis of energy-minimized 



configurations [11] we find that it is necessary to include 
corrections to scaling, but also that the magnitude of the 
corrections are markedly different for different quantities, 
and, furthermore, that the neglect of these corrections is 
a major reason for the differing values for the critical ex- 
ponents in the literature. We find strong evidence for a 
superlinear behavior of yield stress and pressure above 
jamming from the scaling analysis, and also find inde- 
pendent support for this result from pressure data ex- 
trapolated to the limit of vanishing shear rate. We also 
suggest a possible mechanism behind this behavior. 

Following O'Hern et al. [2] we use a simple model of 
bi-disperse frictionless soft disks in two dimensions with 
equal numbers of disks with two different radii in the 
ratio 1.4. Length is measured in units of the diameter 
of the small particles, d s . With r,j the distance be- 
tween the centers of two particles and dij the sum of 
their radii, the interaction between overlapping particles 
is V(rij) = {e/2){\ — rij I dij) 2 . We use Lees-Edwards 
boundary conditions [12] to introduce a time-dependent 
shear strain 7 = fry. With periodic boundary conditions 
on the coordinates Xi and yi in an L x L system, the 
position of particle i in a box with strain 7 is defined as 
Vi = (xi + jyi,yi). We simulate overdamped dynamics 
at zero temperature with the equation of motion [13], 

dt~ C ^ dn +Vllx > 

with e = 1 and C = 1. The unit of time is r = d s /Ce. 
All our simulations at the lower shear rates are from N > 
65536 total particles. 

Our basic scaling assumption describes how different 
quantities, as e.g. shear stress, pressure, potential energy 
and jamming fraction, depend on a change of length scale 
with a scale factor b: 

0(8cb, 7, 1/L) = b-yo^goiScbb^^jb^b/L). (1) 

Here Sqb = <j> — (f>j, yo is the critical exponent of the 
observable O, v is the correlation length exponent, and 
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FIG. 1. Approximate determination of (j>.j and q p from Eq. (2) 
without corrections to scaling. The figure shows pressure ver- 
sus shear rate at several different packing fractions. The pres- 
sure is shown as p/7 ' 3 to make the behavior clearly visible. 
This suggests that p ~ 7 91 with qi = 0.30 at 4>j! = 0.8433. 

z is the dynamic critical exponent. Point J is at Sep = 0, 
7 — > 0, and an infinite system size, 1/L —> 0; the scaling 
relation describes the departure from the critical point in 
these respective directions. 

The above expression may be used as a starting point 
for our analysis. We make use of data obtained at finite 
shear rates and system sizes large enough that finite size 
effects may be neglected — essentially the same approach 
as in Ref. [3]. With b = A f~ 1 / Z in Eq. (1) and qo = Vol zv , 
the scaling relation becomes 

OW,T)~f°9oW/7 1/w ): (2) 

where the scaling function is a function of only a single 
argument. At <fij we have 0(</>j,7) ~ 7 90 which gives 
a simple method for determining qQ and 4>j: Plot O 
versus 7 on a double-log scale for several different <j>. The 
packing fraction for which the data fall on a straight line 
is then our estimated <j>j. Data above and below (pj, 
respectively, should curve in opposite directions. 

We start by applying this simple recipe to the pressure, 
p, and will turn to the shear stress only as the next step. 
Both these quantities are calculated, as in Ref. [2], from 
the elastic forces only. Figure 1 shows pressure versus 
shear rate at several different packing fractions. Antic- 
ipating that the value of q p » 0.3, we plot p/7 3 vs 7 
in order to more clearly differentiate the behaviors near 
<t>j. It is then easy to identify the density with a recti- 
linear behavior, and we find p ~ j qi with q\ = 0.3 at 
4>ji = 0.8433. Data at lower and higher densities curve 
downwards and upwards, respectively. These values 4>j\ 
and q\ are only first estimates of the jamming density 
and the exponent, respectively; our final estimates turn 
out to be just slightly different. 

Figure 2 is the same kind of plot for the shear stress, 
a, and it is immediately clear that these data are not 
directly amenable to the same kind of analysis; there is no 
density with an algebraic behavor across the whole range 
of shear rates. Before presenting our further analyses we 
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FIG. 2. Shear stress a versus shear rate 7 at several different 
densities. Panel (a) shows that there is no density where a 
behaves algebraically across the extended range of shear rates, 
however data across two orders of magnitude of 7 could, to 
a reasonable approximation, be taken as algebraic. In that 
vein, data limited to IO -6 < 7 gives <f> = 0.8416 (crosses) as a 
good candidate for <f>j (cf. Ref. [3]) whereas other ranges of 7 
would give other estimates. From a comparison with p/7 91 = 
const at 4>Ji i n Fig- 1> panel (b) shows the correction term 
a/j qi — g 1 at <f>ji, and it appears that this correction to a 
very good approximation is ~ 7" , which has the same form 
as standard corrections to scaling in critical phenomena. 

note that this provides an explanation for the differing 
values of both jamming density and exponents in the 
literature. In Ref. [3] the jamming density was found to 
be w 0.8415 and the figure shows that data in the range 
10~ 6 < 7 < 10~ 4 would suggest 4> = 0.8416 (crosses) as 
a good candidate for <pj. However, it is clear that data 
at the same density and lower shear rates deviate from 
the algebraic behavior. Similarly, with access to a down 
to 7 = 10~ 7 , cf) = 0.8424 (open circles) would appear as a 
good candidate for <f>j, whereas data in the range 10~ 8 < 
7 < 10~ 6 would suggest 4>j = 0.8433 (solid dots). The 
value of the effective exponent q a also changes: for these 
three different ranges of shear rates we find q a = 0.44, 
0.41, and 0.33, respectively. Note that this explanation 
is at variance with Ref. [8] where the differing exponents 
are attributed to using data from a too large range in <\>. 
That explanation is not applicable here since our analyses 
only consider data right at the presumed <pj. 

As a step towards the final analysis we now consider 
the shear stress at 4>j\ and focus on the deviation from 
the algebraic behavior ~ 7 91 . From Fig. 2(a) we note 
that a/j qi in the limit of low 7 appears to saturate at a 
finite value, 0.005 < gi < 0.006 and so we plot <r/j qi —g\ 
in Fig. 2(b). It is then possible to adjust g\ such that the 
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remainder is algebraic in 7, 

a(0j,o,7)/7 91 =9i+i S ' h , 



(3) 



with the exponent Coq ~ 0.25. 

The importance of this observation lies in the fact that 
standard corrections to scaling modify Eq. (1) to give 
precisely this form [14], 



yo/v _ 



where ho is another scaling function and u) is the correc- 
tion to scaling exponent. Using b — 7~ 1 / z in the above 
then gives 

0(^, 7 )/7 go =go(6(b/j 1/z n+r /z ho(5(b/j 1/z n- (4) 

Equation (3) is just the special case when 5(j) = 0. 

The above analysis of a relied on <\>j\ and q\ deter- 
mined from the pressure data without corrections to scal- 
ing. We now set out to analyze both pressure and shear 
stress directly from the scaling relation, Eq. (4) , that in- 
cludes the correction term, and determine the 4>j, qo, 
1/zv, and lu/z that allow for the best fit to Eq. (4). Here 
go and ho are scaling functions which we approximate 
with fifth-order polynomials in x = 84>/^ 1 l zl ' . The actual 
fits are done by minimizing x 2 /dof with a Levenberg- 
Marquardt method. The number of points in the fits 
range from about 100 to 250 depending on the range of 
data used in the fits. 

In this kind of involved analysis it is crucial to validate 
the results and to that end we use several different crite- 
ria: (i) The first is to check the quality of the fits: Are 
the deviations of the data from the scaling function con- 
sistent with the statistical uncertainties? We use x 2 /dof , 
which should be close to unity to get a quantitative mea- 
sure, (ii) A good quality of the fit does however not by 
itself guarantee that the results are reliable. The sec- 
ond check is therefore whether the fitting parameters are 
reasonably independent of the precise range of the data 
included in the fit. We do this by systematically varying 
both the range of shear rates and the range of densi- 
ties; fixing X — (4> — 0.8434)/7° 26 we use the criterion 
\X\ < A max with X max = 0.2, 0.3, and 0.4. This restric- 
tion does not reflect the size of the critical region but 
rather that the polynomial approximation of the scaling 
function breaks down for too large X. (iii) A final check 
is whether the critical parameters from analyses of dif- 
ferent quantities (here p and a) agree with one another. 

Figures 3 show x 2 /dof and the key fitting parameters, 
ip j, 1/zv, q p and q a plotted against 7 max . For each quan- 
tity the left and right panels are from analyses of pressure 
and shear stress, respectively. First considering x 2 /dof 
in the first pair of panels, we note that the fits are only 
good when the data are taken from a rather restrictive 
interval in <j> around <j>j, \X\ < 0.3. For pressure there is a 
good fit to the data over a very large interval — more than 
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FIG. 3. Results from scaling analyses that include corrections 
to scaling. The left and right panels are from analyses of pres- 
sure and shear stress, respectively. The first pair of panels, 
which give x 2 /Aol, suggest that the analyses are only reliable 
when data is used in a rather restrictive interval of cj> — <j>j, 
\X\ < 0.3. For shear stress one also has to be restrictive in us- 
ing data with larger 7. From the following panels we read off 
(j>j = 0.84347, 1/zv = 0.26, and q p = q a = 0.28. Combining 
the last two (y = qzv) gives y v — y„ — 1.08. 



four decades in 7. For the shear stress the highest shear 
rates should not be used, and reliable results are obtained 
by restricting 7 to 7 < 5 x 10 -5 when X max = 0.2 and 
7 < 1 x 10" 5 for X max = 0.3. 

The next two panels show <pj from pressure and shear 
stress, respectively, in good agreement with one an- 
other; we estimate <f>j = 0.84347 ± 0.00020 in agreement 
with other recent determinations of 4>j from quasistatic 
simulations [11, 15]. Here and throughout, the error bars 
in the figures are one standard deviation whereas the nu- 
merical values give a min-max interval (± three standard 
deviations) for the estimated quantities. To correctly in- 
terpret these figures one should note that the fitted val- 
ues for different 7 max and A max are based on different 
subsets of the same data, and therefore are expected to 
be strongly correlated. The main point is here to check 
how robust the fitting parameters are to changes in the 
precise data set, and the absence of clear trends in the 
results is therefore an encouraging sign. 

We further find l/zv = 0.26 ±0.02 and q = 0.28 ±0.02. 
Combining the two exponents we find y = qzv = 1.08 ± 
0.03 (a strong correlation between q and 1/zv is re- 
sponsible for the error estimate for y). Since y is just 
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FIG. 4. Illustration of results of the scaling analysis. The 
dashed lines are the behaviors at <j>j for p and a, respectively: 
h,i)h qo = go{0) + j uj/z h o (0). 



slightly above unity we have also reanalyzed the pres- 
sure data with the assumption y„ = 1. The fits then 
become considerably worse and we conclude that the 
data is strongly in favor of y p > 1. A similar analysis 
of the shear stress is not conclusive. Using v = 1.09 
from Ref. [11] the dynamic critical exponent becomes 
z = 3.5 ± 0.4. The correction to scaling exponent (not 
shown) is oj/z = 0.29 ± 0.03, or ujv = 1.10 ± 0.06, which, 
again using v = 1.09 [11], gives cj = 1.0 ± 0.1 in good 
agreement with Ref. [11]. 

The analyses of both pressure and shear stress work 
nicely when corrections to scaling are included. A draw- 
back with including the corrections is — beside the more 
difficult analyses — that it is no longer possible to deter- 
mine <j)j directly from a simple plot as in Fig. 1. The 
most direct way to illustrate the determination of <f>j is 
shown in Fig. 4 which displays p/ A f Qp and a/j q " against 
7^/ z , now with linear scales on both axes. Data at <j)j 
should then fall on a straight line. Note the very different 
size of the corrections, given by the slopes of the data. 

For well above <pj the pressure decays algebraically in 
7 and this gives a means to determine the limiting value 
p((f>, 7 — > 0). If we can get reliable values, p(4>, 7 — > 0), at 
densities sufficiently close above it should be possible 
to get another determination of y p , independent of the 
scaling analysis. Fig. 5 shows some of our finite-7 data 
together with such extrapolated values for densities down 
to <j> = 0.848. Fitting to the five points from 4> = 0.848 
through 0.856 (0.5% through 1.5% above <j)j) we find y = 
1.09 ±0.04 shown by the solid line, in excellent agreement 
with y = 1.08 from the scaling analysis. (The inset of 
Fig. 5 shows how y depends on the assumed <j)j.) Similar 
results, y v k, 1.1 have also been found before [16, 17]. 

The above results point to a good agreement between 
the exponent obtained from the scaling analyses on the 
one hand, and the 7—^0 limit of the pressure above 
<pj on the other. This is entirely in accordance with ex- 
pectations from critical scaling. This is in contrast to 
the claim in Ref. [8] that the critical region is extremely 
narrow and doesn't include densities away from <j)j in the 
limit the yield stress is there taken to be governed 

by a different regime with a different exponent, y a = 3/2. 

With the result y « 1.1 from two different analy- 
ses it becomes important to try and reconcile this with 
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FIG. 5. Alternative determination of the exponent y p . The 
open circles are p(</>,7 — > 0) from extrapolations of p(<f>,-y). 
Assuming (j>j — 0.84347 the exponent becomes y — 1.09, 
shown by the solid line. The dashed line corresponds to y = 1. 
The inset shows how y depends on the assumed <f>j. 



the well established linear increase of the pressure when 
marginally jammed packings are compressed above their 
respective jamming densities [2, 9]. We speculate that the 
reason for this is that the ensemble of configurations de- 
pends in a non-trivial way on (f> in the vicinity of <f>j, and 
that this is so since the dynamic process that generates 
this ensemble is itself very sensitive to <fi. It is then rele- 
vant to consider the behavior in the quasistatic limit and 
to recall that the average time needed for the minimiza- 
tion of energy in quasistatic simulations diverges as <pj 
is approached from above or from below. (This parallels 
the more rapid jumping between jammed and unjammed 
states reported in Ref. [15].) A dramatic change of the 
dynamical process that generates the ensemble suggests 
that the ensemble itself would depend on ^ in a non- 
trivial way. 

To conclude, we have shown that pressure and shear 
stress from shearing simulations are entirely consistent 
with the assumption of a critical behavior when correc- 
tions to scaling are included in the analysis. We find 
4>j = 0.84347 ± 0.00020 and that at <j>j, both p and a 
scale as 7 9 with q = 0.28 ± 0.02. In the limit 7^0 both 
p and a vanish as (<f> — <j>j) v with y — 1.08 ± 0.03. 
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